#include "SPDE_One_Mode.h"

// --------------------------------------------------------------
SPDE_One_Mode::SPDE_One_Mode(int ip,int im, int in, int L, double idt,double iXYmax):
p(ip), m(im),n(in),LatticeSize(L),XYmax(iXYmax)
{
  // BASIC BLOCK ------------------------------
  // Input#0 = T (K)
  // Input#1 = W (nm/s)
  // Input#2 = A (nm/s)
  // Input#3 = Freq
  int idim[4] = {1,1,1,1};
  // Output#0 = z_p (nm)
  // Output#1 = (z_p)^2 (nm^2)
  int odim[2] = {1,1};
  block::initialize(4,idim,2,odim,idt);

  // STATES -----------------------------------
  z = block::outputs[0];
  z_next = new double[1];

  // MODEL PARAMETERS -------------------------
  lambda = 0.0;

  // NUMERICAL SOLUTION PARAMETERS ------------
  // For Analytical Expression
  sqrt_dt = sqrt(block::dt);

  // Setup random number generator
  gsl_rng_env_setup();
  noise_type = gsl_rng_default;
  noise = gsl_rng_alloc(noise_type);
  gsl_rng_set(noise,time_seed());
  
  // INITIALIZATION -------------------------------
  reset();
}
// --------------------------------------------------------------
SPDE_One_Mode::~SPDE_One_Mode(){
  delete []z_next;
 gsl_rng_free(noise);
}
// --------------------------------------------------------------
void SPDE_One_Mode::reset(){
  block::time = 0.0;
  // Reset states
  z[0] = 0.0;
  z_next[0] = 0.0;
  // Reset outputs
  outputs[1][0] = getR2();
}
// --------------------------------------------------------------
void SPDE_One_Mode::setPara(int para_id,int para_val){
  switch(para_id){
    case 0:
      p = para_val;
      break;
    case 1:
      m = para_val;
      break;
    case 2:
      n = para_val;
      break;
  }
}
// --------------------------------------------------------------
void SPDE_One_Mode::update(double sysTime){
  if(sysTime<time){
    cerr << "Error@SPDE_One_Mode::update(), sysTime = " << sysTime << ", time = " << time << "\n";
    exit(1);
  }

  // Update states  ====================================
  while(time+dt<=sysTime || time == 0){
    if(sysTime!=0){
      z[0] = z_next[0];
    }
    // Calculate x_{k+1} ---------------------------
    T      = inputs[0][0];
    W      = inputs[1][0];
    A      = inputs[2][0];
    Freq   = inputs[3][0];

    if(T>1500 || T<273){
      cerr << "Invalid value of T: " << T << "\n";
      exit(1);
    }

    if(W>10.0 || W<0.0){
      cerr << "Invalid value of W: " << W << "\n";
      exit(1);
    }

    if(A>W && A<0.0){
      cerr << "Invalid value of A: " << A << "\n";
      exit(1);
    }

    if(Freq<0){
      cerr << "Invalid value of Freq: " << Freq << "\n";
      exit(1);
    }
      
    // UPDATE MODEL PARAMETERS ---------------------
    double c2 = eq_c2(T,W);;
    double sigma2 = eq_sigma2(T,W);
    lambda = -4*c2*(m*m+n*n)*pi*pi/(XYmax*XYmax);

    // UPDATE STATES -------------------------------
    // Analytical expression ==============================
    double temp1 = exp(lambda*dt);
    double temp2 = sqrt(sigma2*(exp(2.0*lambda*dt)-1.0)/(2.0*lambda));
    double temp3 = (temp1-1.0)/lambda;
  	
    switch(p){
      case 1:
	if(m>0&&n>0){
	  z_next[0] = z[0]*temp1+temp2*gsl_ran_gaussian_ziggurat(noise,1.0);
	}
	break;
      case 2:
	if(m==0&&n==0){
	  // double w200 = W*XYmax+A*XYmax*(1-cos(2*Freq*pi))/(2*Freq*pi);
	  double w200 = W*XYmax;
	  z_next[0] = z[0]+w200*dt+gsl_ran_gaussian_ziggurat(noise,1.0)*sqrt(sigma2)*sqrt_dt;
	}
	else{
	  if(n==0 && Freq!=m){
	    // double w2m0 = SQRT_TWO*A*XYmax*Freq*(cos(2*Freq*pi)-1)/(2*pi*(m*m-Freq*Freq));
	    double w2m0 = 0.0;
	    z_next[0] = z[0]*temp1+
	      temp2*gsl_ran_gaussian_ziggurat(noise,1.0)+temp3*w2m0;
	  }
	  else{
	    z_next[0] = z[0]*temp1+temp2*gsl_ran_gaussian_ziggurat(noise,1.0);
	  }
	}
	break;
      case 3:
	if(m>0){
	  if(n!=0){
	    z_next[0] = z[0]*temp1+temp2*gsl_ran_gaussian_ziggurat(noise,1.0);
	  }
	  else{
	    if(m==Freq){
	      double w3k0 = SQRT_TWO*A*XYmax/2;
	      z_next[0] = z[0]*temp1+temp2*gsl_ran_gaussian_ziggurat(noise,1.0)+temp3*w3k0;
	    }
	    else{
	      // double w3m0 = SQRT_TWO*A*XYmax*m*sin(2*Freq*pi)/(2*pi*(Freq*Freq-m*m));
	      double w3m0 = 0.0;
	      z_next[0] = z[0]*temp1+
		temp2*gsl_ran_gaussian_ziggurat(noise,1.0)+temp3*w3m0;
	    }
	  }
	}
	else{
	  z_next[0] = 0.0;
	}
	break;
      case 4:
	if(n>0){
	  z_next[0] = z[0]*temp1+temp2*gsl_ran_gaussian_ziggurat(noise,1.0);
	}
	else{
	  z_next[0] = 0.0;
	}
	break;
    }

    // ==========================================
    if(time==0)
      break;
    else{
      time=time+dt;
    }
    // ==========================================
  }// end of while

  // OUTPUTS -----------------------------------
  outputs[1][0] = getR2();
  return;
}
// --------------------------------------------------------------
double SPDE_One_Mode::getR2(){
  return (z[0]*z[0]);
}
// --------------------------------------------------------------
double SPDE_One_Mode::eq_c2(double T,double W){
  // INPUTS:
  // T: k
  // W: nm/s
  //
  // OUTPUTS:
  // C2: nm^2/s
  //
  //  double pc1 = 7.456e-8;
  //  double pc2 = -4.039e-7;
  //  double pc3 = 8.092e-7;
  //  double pc4 = -7.376e-7;
  //  double pc5 = 7.579e-7;
  //  return pc1*W*W*W*W+pc2*W*W*W+pc3*W*W+pc4*W+pc5;
  //  return k_c*exp(k_B*T/E_c)*pow(W,-a_c)/(LatticeSize*LatticeSize);
  double pc1 = -10.53;
  double pc2 = 12.23;
  double pc3 = -4.953;
  double pc4 = 0.8804;
  return pc1*W*W*W+pc2*W*W+pc3*W+pc4;
}

double SPDE_One_Mode::eq_sigma2(double T,double W){
  //  double ps1 = 1.224e-4;
  //  double ps2 = 2.871e-6;
  //  return ps1*W+ps2;

  //  double index = (a_t+k_t*W)*T-a_v-k_v*W;
  //  return pi*pi/(LatticeSize*LatticeSize)*(1+exp(index));

  // 
  double ps1 = -15.98;
  double ps2 = 16.39;
  double ps3 = -5.634;
  double ps4 = 0.8773;
  double ps5 = 0.01836;

  return ps1*W*W*W*W+ps2*W*W*W+ps3*W*W+ps4*W+ps5;
}
// -------------------------------------------------------------
